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ABSTRACT 



The M-Layer program tracks the constant phase lines Im{D(q)} = 0 and looks for 
their intersections with the lines Re{D(q)} = 0 for the locations of the zeros of the mode 
function D(q). These two types of constant phase lines are tracked and plotted over a 
search region which contains modes having a range attenuation rate of no more than 5 dB 
per km. Several new parameters for use in mode search are deduced from the results and 
some old ones are verified. Future studies in waveguide mode propagation theory 
pertaining to atmospheric ducts may benefit from this work. An improved mode search 
strategy is also proposed. 
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L INTRODUCTION 



The M-Layer program developed by NOSC (Naval Ocean System Center) was 
documented by Yeoh [Ref. 1] at NFS. It was later extensively revised by Lee and Han 
[Ref. 2] for improved accuracy and speed. The new FORTRAN code is now identified 
as the NFS version [Ref. 3] under the auspices of NRaD (Naval Command, Control and 
Ocean Surveillance Center, RDT&E Division), the current name for NOSC. In this 
version, the logical structures and numerical algorithms for following the constant phase 
lines of the mode function and for computing the mode locations were almost completely 
rewritten. The overall plan to first partition the search region into rectangular areas and 
then circulate around these rectangles to search for the modes, i.e., the mode search 
protocol, was left unchanged. Lee and Han [Ref. 2] recommended that a more direct 
mode search protocol should be devised to locate the modes more expediently. If possible, 
such an approach should locate the modes in the order of ascending range attenuation 
rates. 

To design a more efficient mode search protocol, a better understanding of the mode 
function beyond the known locations of the modes is required. In this thesis, the 
analytical property of the mode function is investigated. Frograms are written to track and 
plot the constant phase lines along which the mode function is either real or purely 
imaginary. From the results, a new strategy to locate the modes is recommended. 
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In the remainder of this chapter, the background theory and some of the notations 
used in this thesis are introduced. The basic theory presented in Sections A and B follows 
Ref. 3 closely. The results of this research are presented in Chapter II. In Chapter III, a 
new mode search protocol is proposed based on the findings of this thesis. The relevant 
parameters for use in this new approach are also discussed. 

A. THE WAVEGUIDE MODE THEORY OF PROPAGATION 

Trapping of electromagnetic (EM) waves in the modes supported by a duct is the 
dominating factor in over-the-horizon propagation. The computer program M-Layer 
searches for these modes and computes the electric field from the Hertz vector. Assume 
a vertical electric dipole 



J = 47tzfi(x)5(y)«(z-zp 

at a height z=Zj above the surface of the earth. Following Freehafer [Ref. 4], the Hertz 
vector of the EM fields can be written, in the cylindrical coordinates (p, <J>, z), as a sum 
of contributions from individual modes [Ref. 1]: 

n(p,z) = -niY, (1) 

m 

where is the wavenumber of the m-th mode and is independent of the coordinates; 
Zj is the height of the transmitter above the ground, Hq^^^ is the Hankel function of the 
second kind, which represents an outgoing wave in the radial direction when the time 
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dependence is adopted, and is the height-gain function of the m-th mode, 
normalized so that the integral of its square over all heights equals unity when either an 
electric or a magnetic dipole source is used. The height-gain function satisfies the 
equation 



where k is the free- space wavenumber, m(z) is the modified index of refraction, and 
m^(z) is approximated with a continuous, piecewise linear profile having I layers: 

m\z)=m^+a^(z-z^) z^^z^z^^^, l^i^L 

In practice, m(z) deviates only slightly from unity in the troposphere. The modified 
refractivity M(z) = [m(z)-l]xlO^, as shown in Fig. 1, is commonly used. The slope of 
M(z) is approximately a -x 10^/2. 

In the i-th layer, the height-gain function is given by 

where q^^ is a dimensionless linear function of height z with k, m^, and nr- as 
parameters: 
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The functions k^(qj^) and function Ai: 



*i(0=-A12)«‘^i(-9^e/^’''’) 



( 6 ) 



and 






(7) 
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For z > ^i^iy rn^(z) is extended continuously upward at a constant slope 
conimensurate with the effective radius of the earth. This slope equals 2.36x10”^ (m~b 
for the four-thirds effective earth radius model. Thus, the height-gain function g^^ is again 
given by Eqs. (4) through (7) for z > Zj+j, with the constant set to to 

represent an upward going wave as z becomes large. Beneath the "flattened" earth surface 
2 

where z < z^ =0, m (z) is deduced from the dielectric constant and the conductivity of 
sea water, which are set to 80.8869 and 4.64 (Si/m) respectively. Here the Hertz vector 
is represented by a downward propagating plane wave in this region. 

The conditions that the tangential components of the electric and magnetic fields 
are continuous across the layer boundaries determine the coefficients C-. Starting from 
the top layer down, or "integration-down" [Ref. 1], every C- is uniquely determined by 
at z = z^^j. Thus, a set of all Cj coefficients is determined without considering the 
boundary conditions at z = z^. On the other hand, starting from the lowest boundary at 
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z = Zj up to z = Zj, or "integration-up'', a second set of C- coefficients is determined 
without applying the boundary conditions at z = Zj^j. That these two sets of C- 
coefficients are identical is a manifestation that is the wavenumber of a mode. 
Mathematically, this consistency condition, sometimes called the guidance condition 
[Ref. 5], can be expressed as the vanishing of a determinant D called the mode function 
[Ref. 6]. The elements of this determinant consist of linear combinations of k^Cq^ni) and 
k 2 (qmi) and their derivatives with respect to height at the boundaries. The M-Layer 
program searches for the zeros of the mode function to find the wavenumbers. Once a 
wavenumber is obtained, the height-gain function of the corresponding mode can be 
computed. Note that the normalization condition on gj^ and the boundary conditions, 
together with the C* coefficients, determine the B- coefficients to within an overall sign. 
This sign need not be resolved because only products in the form of 
each mode contribute to the Hertz vector. 

Since the mode function D depends on the wavenumber p^^ only through the values 
of q^^ at the layer boundaries, it is more convenient to consider D as a function of the 
variable q given by 



and to search for the zeros of D(q) in the complex q plane. The m-th zero is designated 
as qj^ and is called a q-eigenvalue and Pj^ is then deduced from q^^ by inverting Eq. (8) 




( 8 ) 
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for p . In the NPS version, is ordered in ascending attenuation rate of the mode, which 
is approximately proportional to the imaginary part of P^^. 



B. MODE SEARCH PRINCIPLE 

M-Layer searches for modes which have attenuation rates below a specified value. 
At ground ranges more than several wavelengths away, this attenuation rate is 
proportional to the imaginary part of the wavenumber Pj^ as can be seen from Eq. (1). 
In the upper complex q plane, for values of q such that 




« 1 , 



( 9 ) 



mj-p/k is approximately proportional to q because the absolute value of m^ and, hence, 
that of p/k, are both nearly unity. Thus, the limit on attenuation rate imposed on the 
imaginary part of p can be translated approximately into an upper bound ^ on the 
imaginary part of q which defines a strip in the complex q plane called the search region. 
This region is covered with layers of identical square meshes whose sides are parallel to 

9 

the imaginary q axis and have a length equivalent to an attenuation rate of 1/32 dB per 



^ Since the absorption is small in air, the modified index of refraction m(z), and 
hence a-, are considered as real quantities in the following discussions. In the actual 
FORTRAN code, they are declared as complex variables. 

9 

This is the default value of the NPS version which can be adjusted through editing 
the variable "dmesh" in an ASCII input file. 
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kilometer. The lower edges of the meshes in the lowest layer sit along the real q axis , 
The meshes in the layer just below the top one contain the upper bound of the search 
region, with the top layer providing some allowance for numerical inaccuracy. The 
program tests each mesh square for the presence of zeros of the mode function D(q). 

To search through all the meshes, the program first divides the mesh covering the 
search region into ’’contour rectangles” with equally spaced vertical lines parallel to the 
imaginary q axis. These vertical lines contain the edges of stacked square meshes and are 
separated by a distance 160 times the side of a mesh. The search commences at the top 
left comer and moves counterclockwise around each ’’contour rectangle” and begins with 
the one whose left edge is defined by 




=2xl0-‘d^. 



( 10 ) 



where is the minimum value of the modified index of refraction profile and 
is shown in Fig. 1. The justification for this choice to begin the mode search is as 
follows: From Eq. (5) and optics, for a wave to propagate freely in space, the dispersion 
relation of the Maxwell equations requires that k^m^(z) > assuming that both 
quantities are real. Thus, the smallest P ^ to support a trapped wave will occur when 



'i 

The NOSC version extends the lower edge slightly below the real q-axis. This 
causes problems in some situations [Ref. 2]. 
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^ m exceeds the minimum of nr(z) . It is not surprising that an argument based 
on optics works within a waveguide mode formulation which is low frequency in 
principle. Lee [Ref. 7] has demonstrated that the earth-flattening approximation actually 
links Mie’s low frequency oriented spherical harmonics to Fock’s high frequency 
diffraction theory through the uniform asymptotics. The mere introduction of an unlimited 
ground range into the Maxwell differential equations incorporates the global feature of 
the radius of the earth into the local equations. 

After the search over the initial rectangle is completed, the program goes on to 
search the neighboring rectangle to the left. If a specified number^ of consecutive 
rectangles of decreasing real q values have been searched without turning up any zero of 
D(q), the program changes direction and starts to search the rectangles to the right of the 
initial "contour rectangle" one by one, with increasing real part of q. After failing 
consecutively to locate any q-eigenvalue again over the specified number of rectangles, 
the program assumes that no more zeros of D(q) exists in the search region. The mode 
search is considered complete and the procedure is terminated. If the array for storing the 
q-eigenvalues is filled up before the search is completed, the search is terminated with 
an error message. 



^ Note that assuming the same real part, an imaginary component of 0^ reduces the 
real part of 0"^^ which may enable the wave to propagate in the z-direction. 

^ This number of consecutive rectangles is an adjustable input variable "nstop" in the 
NFS version. The default value is 2. 
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The search for zeros of D(q) makes use of the fact that a real valued function 
changes sign when it crosses a simple zero. Since a zero of the complex valued function 
D(q) is where both its real part and imaginary part vanish, a necessary condition for a 
point to be a zero is that it is the intersection of two curves defined by Im{D(q)} = 0 
and Re{D(q)} = 0. M-Layer moves along each side of a "contour rectangle" while 
searching for a sign change in Im{D(q)} across an edge of a mesh bordering the side of 
this rectangle to determine that a line of Im{D(q)} = 0 has been encountered. The search 
then follows this line into the meshes within the "contour rectangle", checking each mesh 
to see if a curve Re{D(q)} = 0 enters the mesh being inspected. All these steps make use 
of the assumption that the zeros of D(q) are simple. Once both the curve Im{D(q)} = 0 
and the curve Re{D(q)} = 0 are found to be present within a mesh, the locations of their 
possible intersections are estimated. 

The rule for sign change becomes inconclusive if a zero of Im{D(q)} = 0 or if 
Re{D(q)} = 0 happens to fall on a comer of a mesh, and a remedy is required. In the 
NFS version, whenever the real part or the imaginary part of D(q) vanishes on a comer 

_C9 

of a mesh, the phase angle of D(q) is rotated by 2 radians. This maneuver effectively 
shifts q by a small amount to resolve the sign ambiguity. 

To estimate the locations of zeros of D(q) within a mesh, D(q) is assumed to be 
well approximated in this mesh by its four-term Taylor series expansion. The unshifted 
values of D(q) at the mesh comers determine this cubic polynomial uniquely. Cardan’s 
formulas [Ref. 8] are used to locate the zeros of this polynomial before retaining only 
those lying within the mesh. 
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n. MODE FUNCTION IN THE CXDMPLEX PLANE 



A. MODE LOCATIONS 

M-Layer searches for the zeros of the mode function D(q) in the complex q plane 
by tracking the constant phase lines of D(q) along which the mode function is either real 
or purely imaginary. The zeros found are called the q-eigenvalues as defined by Eq. (8) 
and are denoted as qj^. These eigenvalues are saved in an ASCII file and are utilized later 
for height-gain function computations. In order to design a strategy for locating these 
zeros, the mode locations are plotted. 

It is clear from Eq. (8) that the variable q varies with a the slope of m (z) in the 
lowest (first) layer. Since depends on how the continuous, piecewise linear 
approximation to m (z) is made, while both m^ and are, in principle, dependent only 
on the actual profile, q will vary strongly with a ^ and fail to provide information on the 
wavenumber of the mode directly. A more suitable variable to use is, from Eq. (8): 




( 11 ) 



Since q is given the variable name of q^j and is given the variable name of t^ 

in the M-Layer FORTRAN code, this variable is identified as q^Aj throughout this 
thesis. Since m^ is real and both it and p/k deviate from unity only slightly for all cases 
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investigated, is approximately 2(mj - ^/k). Thus, its imaginary part is directly 

proportional to -P, the attenuation rate of the mode. Furthermore, q^j/tj does not depend 
explicitly on a^. Removing the tj factor from q^ j makes it possible to compare results 
from different evaporation duct profiles meaningfully. 

Plots of the q-eigenvalues as individual points in the complex plane are 

included in Appendix A. For the 20 meter duct, which will be used as the representative 
case for discussions in this thesis, a line linking all the eigenvalues in the order of 
increasing attenuation rate is shown in Fig. 2. Because of the long excursion between 
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many pairs of consecutive modes in the upper part of this figure, it appears that the 
intuitive approach of starting with the mode of the lowest attenuation and continually 
looking for the mode with the next higher attenuation rate may not be the most efficient 
strategy. Furthermore, there seems to be a fork near the middle of the figure which 
suggests that there can be more than one constant phase line passing through all the zeros 
of the mode function. These problems suggest that, in order to design an efficient mode 
search strategy, the behavior of the mode function in the complex ^/t^ plane has to be 
investigated more thoroughly. 

B. PHASE LINE TRACKING 

The M-Layer program follows a line along which Im{D(q)} = 0 until a zero is 
encountered at its intersection with another line along which Re{D(q)} = 0. Since there 
is no way to determine the phase of the constant phase line linking two adjacent zeros a 
priori^ following these two types of somewhat arbitrary but easy to compute phase lines 
which have constant phases of 0 or ii, and tt/ 2 or 3n/2, respectively, are still the best way 
to program a computer to locate the zeros systematically. It is thus necessary to find out 
the actual distribution of these constant phase lines in the complex plane. 

1. Downward Tracking 

In the original design, along the real q axis, M-Layer divides the search region 
into "contour rectangles," each of which spans 160 meshes horizontally. From the 
observed locations of the modes, it is found desirable to track and plot the constant phase 
lines along Im{D(q)}=0 and Re{D(q)} = 0 over a little more than six "contour 
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rectangles", roughly three to the right and three to the left of the starting real q coordinate 
given by Eq. (10), which is a variable called q^^^^ in M-Layer. The subroutine FNDMOD 
and FZEROX are modified to track these constant phase lines. A listing of these modified 
routines for tracking the lines Im{D(q)) = 0 beginning from the top edge and moving 
downwards into the search region is included in Appendix B. Their flow charts are given 
in Appendix C. Only minor changes are required to track the type of lines Re{D(q)} = 0, 
or to track these lines from the bottom of the search region upwards. The program first 
searches for a sign change in Im{D(q)) or in Re{D(q)} along the top edge of the search 
region, starting at the coordinate -512 mesh sizes from until a distance of 1024 
mesh sizes is covered. When a sign change is observed, a constant phase line is 
recognized as passing through the particular mesh square and the program writes the 
complex q value of the lower-left comer of this mesh square into an ASCII file identified 
by whether the mode function is real or imaginary along the line, and the order this line 
is found. The program then moves from the top edge downwards to follow this constant 
phase line until it exits the search region. The q values of the lower-left comer of the 
mesh squares along this phase line are also written into the same ASCII file to be plotted 
later using MATLAB/386. The constant phase lines for the mode functions of the 2, 4, 
6, 8, 10, 20, 30 and 40 meter evaporation ducts are obtained and plotted. They are 
included in Appendix D. Tracking and plotting these constant phase lines are extremely 
time consuming. The CPU time required to track the desired constant phase lines for each 
duct using an Intel 80486 based PC running at a 33 MHz clock rate is listed in Table 1. 
The time required to plot those lines for each duct using an Intel 80386 based PC mnning 
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at a 16 MHz clock rate is also listed. 



TABLE 1 

CPU Time for Tracking and Plotting Constant Phase Lines 



duct 


tracking 


tracking 


plotting 


height 


Im{D(q)) = 0 


Re{D(q)} =0 


both lines 


(m) 


(hr:min:sec) 


(hr:min:sec) 


(hr:min:sec) 


2 


1:47:39 


1:46:23 


0:21:28 


2 


2:47:41 


2:48:50 


0:23:17 


8 


3:43:54 


3:45:11 


0:25:42 


8 


4:25:49 


4:23:41 


0:29:34 


10 


6:08:14 


6:11:03 


0:32:25 


20 


8:26:58 


8:24:22 


0:34:21 


30 


9:08:14 


9:06:45 


0:36:49 


40 


8:39:04 


8:37:54 


0:36:53 


Total 


45:07:33 


45:04:09 


4:01:29 



The constant phase line plot for the mode function of the 20 m evaporation 
duct is given in Fig. 3. The solid lines represent those having Im{D(q)} = 0; the dotted 
lines represent those having Re{D(q)} = 0. Every intersection of these two types of phase 
lines is the location of a q-eigenvalue scaled by t^, thus indicating the presence of a 
mode. Note that every solid line intersects with a dotted line except for many of those 
running from the top through the bottom of the search region. None of the solid lines 
intersect with other solid lines, neither do the dotted lines. This is seen clearly in Fig. 4, 
which magnifies the lower center portion of Fig. 3. One feature that can be recognized 
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immediately in all the plots in Appendix D is that the starting real q coordinate given by 
Eq. (10), which is marked by an "x" along the top edge of each plot, correlates very 
closely with the mode of minimum attenuation. 

2. Upward Tracking 

One mode of low attenuation rate is missing from Fig. 2 compared to those 
found in Ref. 2. Several of them are also missing from those cases of greater duct heights 
plotted in Appendix D. It is evident that tracking the constant phase lines from the real 
qil/tj axis upwards is necessary. For the 10 meter and the 20 meter ducts. Fig. 5 shows 
several constant phase lines starting out of the lower limit of the search region and 
returning to the lower-half complex q^ j/t^ plane. Fig. 6 shows the presence of such lines 
for the 30 meter and the 40 meter ducts. These constant phase lines have not been 
observed in those cases of lower duct heights. 
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and 20 meter (bottom) ducts. 
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m. ANALYSIS AND CONCLUSIONS 



A. MODE SEARCH PARAMETERS 

The implementation of a mode search procedure requires several parameters. The 
search region has to be defined first. In M-Layer, the desired maximum range attenuation 
rate of the modes which support wave propagation is treated as an input parameter called 
^loss’ conventionally in dB per km. This parameter determines the region over 
which the modes are to be searched as explained in Chapter I. In the complex q plane, 
under the assumption that both m^ and p/k of interest are close to unity, the search region 
is bounded from above by the FORTRAN variable q^^p given by 



^top 



*xlogio« 



®loi« • 



( 12 ) 



The location at which to start the mode search, is another parameter 

determined by M-Layer. Based on optical considerations, the M-deficiency, of 

Fig. 1, which is the amount of decrease of the modified refractivity from its value on the 
sea surface to its minimum value in the air, is a measure of the capability of the duct to 
trap electro-magnetic waves. As explained also in Chapter I, this quantity determines the 
location for the program to start searching for modes. The validity of the optical argument 
and the usefulness of Eq. (10) have been proved by this work, as pointed out in Chapter 
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n. 



To search for zeros of the mode function, M-Layer first covers the search region 
with identical mesh squares. It then follows the lines of the type Im{D(q)) = 0 through 
each mesh square, checking for indications that a curve Re{D(q)} = 0 is entering the 
same mesh so that an intersection is possible. The term "mesh size" will denote the size 
of the edges of the mesh squares. It is the size of the step taken by the program to 
advance through the search region. It also determines the initial resolution of the location 
of a mode. The choice of the mesh size should strike a balance between the desire for a 
speedy completion of the search process and the requirement in mode locating accuracy. 
As reported in Ref. 2, for all the evaporation ducts considered, the choice of a mesh size 
equivalent to an attenuation rate of 1/32 dB per km appears to be optimal, that is, all 
modes can be found for all cases investigated in Ref. 2 without experiencing 
extraordinarily long computation time. Consider this mesh size as the default and call it 
q^ in the complex q plane. Under the same assumption for deducing Eq. (12), the value 
of q^/tj at 9.6 GHz equals 3.58x10”^. Figure 7 shows the separation between two 
adjacent constant phase lines of the type Im{D(q)} = 0, indicated in the figure as solid 
lines, along the top edge of the search region for the 20 meter duct. The minimum of 
these separations in terms of is about 2.8x10 , which corresponds to a spacing 

of a little less than eight mesh squares apart. The minimum separation between adjacent 
Im{D(q)} = 0 and Re{D(q)} = 0 lines is thus about four meshes. Constant phase line 
separations for all evaporation ducts considered are included in Appendix E. Note that the 
minimum separation is almost constant for ducts higher than 20 meters. It increases 
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(20 m duct). 
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slightly as duct height is decreased. 

As the program follows an Im{D(q)) = 0 line into the search region, a neighboring 
Re{D(q)} = 0 line moves closer and enters into the same mesh square. This causes the 
program to invoke the root finding routine to determine the possible locations of zeros 
of the mode function within this mesh square. If the mesh size is too large and the 
Re{D(q)} = 0 line enters the mesh square long before the actual intersection takes place, 
the root finding routine will be prematurely invoked many times and the probability of 
producing false modes is increased. This explains why a larger mesh size sometimes will 
increase the execution time. 

Given these three parameters q^^p, and q^, a mode search strategy which does 
without the "contour rectangles" is proposed. It should improve the efficiency of M-Layer. 

B. MODE SEARCH STRATEGY 

From the constant phase line plots of Appendix D, a strategy to search for the mode 
can be drawn: move along the top edge starting at Search first toward either the left 
or the right, then reverse course to search along the other direction. After the search along 
the top edge is completed, search the lower edge from one end to the other. In what 
follows, the implementation of this strategy will be discussed. 

1. Track Termination 

The tracking of an Im{D(q)) = 0 line naturally ends if the top edge or the 
bottom edge of the search region is reached. On the other hand, the search region is 
unbounded to the right and left of It is convenient to retain the feature in the 
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original program to set a limit on the number of steps allowed to follow a constant phase 
line. From Fig. 3, this limit can be set to 2.5 times the number of mesh squares between 
the vertical limits of the search region. This number equals 2.5x32xa|Qgg. 

2. Search Termination 

M'Layer has to determine that no more modes within the search region is to 
be found and to terminate the search for modes. When searching along the top edge of 
the search region for the constant phase lines Im{D(q)} = 0, the separation between the 
real parts of the mode eigenvalues is of interest Figure 8 shows the separation in the real 
part of neighboring q-eigenvalues scaled by tp plotted against their locations q^/tj along 
the real axis for the case of the 20 m duct Similar plots for all duct heights are 

grouped in Appendix F. Excluding the lowest point which involves the mode obtained via 
searching along the lower edge, the separation between neighboring modes is erratic with 
an upward trend away from the center of the figure, which is close to the search starting 
position. 

The increase in distance between two modes towards the end of the range 
searched makes it difficult to implement an adaptive mechanism to terminate the search. 
On the other hand. Fig. 3 shows that almost all phase lines entering the search region 
from the top that contain a mode within this region are bunched together. Therefore, the 
parameter set equal to four and may be adjusted, can be used to stop the search 
along the top edge when four consecutive constant phase lines tracked turn up no mode. 

The search along the real q axis can be confined to within the end points of 
the search along the top edge. 
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Figure 8. DiJferencc between consecutive Re(q^tj) values (20 m duct). 
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3. Track Duplication Avoidance 

When the program searches step by step along the top or bottom edges for 
sign changes in Im{D(q)} to start tracking the constant phase line, the exit of a constant 
phase line from a previous track will be encountered. Entering into the search region at 
such a location will simply retrace a constant phase line which has already been examined 
for the existence of a mode. Thus, the exit locations of the constant phase lines must be 
recorded and checked whenever a sign change in Im{D(q)} is found before deciding to 
follow the constant phase line into the search region. Along the top edge, a record of four 
updated most recent exit locations from the top edge should be adequate. The locations 
of exits from the bottom should be kept for use during the search along the lower edge 
of the search region. For this record, a dimension of 256 should be adequate for the 
current cases. But this dimension should be adjusted if other types of ducts are studied. 
A record of four of the most recent exit locations out of the lower edge should also be 
kept and checked to avoid duplicate tracking of the constant phase lines. 
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APPENDIX A: MODE LOCATIONS IN THE COMPLEX PLANE 
This appendix contains figures of mode locations of evaporation ducts of 2, 4, 6, 
8, 10, 20, 30, and 40 m heights plotted in the complex q^j/t^ plane for easy comparison. 
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Hgure A.l Qu/tj mode locations (2 m duct). 
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Rgure A.2 mode locations (4 m duct). 
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Figure A.3 jA j mode locations (6 m duct). 
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Bgure A.4 Qu/tj mode locations (8 m duct). 
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Figure A.5 QiiAj mode locations (10 m duct). 
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Figure A.6 quAj mode locations (12 m duct). 
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Figure A.7 qu/t| mode locations (14 m duct). 
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Figure A.8 mode locations (16 m duct). 
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Figure A.9 qijAj mode locations (18 m duct). 



37 




Figure A.IO qiiAj mode locations (20 m duct). 
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Figure A.ll qjj/tj mode locations (22 m duct). 
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Figure A.12 j/tj mode locations (24 m duct). 
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Figure A. 13 qji/tj mode locations (26 m duct). 
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Figure A.14 mode locations (28 m duct). 
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Figure A.15 QiiAj nicxle locations (30 m duct). 
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Figure A.16 q||/t| mcxic locations (32 m duct). 
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Figure A.17 Qu/ti mode locations (34 m duct). 
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Figure A.18 qij/tj mode locations (36 m duct). 
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Figure A. 19 QuAj mode locations (38 m duct). 
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Figure A.20 qjiAi mode locations (40 m duct). 
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APPENDIX B: SUBROUTINE FNDMOD AND FZEROX 



This appendix contains the listing of the subroutines FNDMOD and FZEROX, 
modified from the NPS version of the M-Layer FORTRAN code to track the constant 
phase-lines and to locate the modes. 
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Line# 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 



Source Line Microsoft FORTRAN Optimizing Compiler Version 5.00 
subroutine fndmod(aloss,dmmin,tl ,dmesh,filemmrmode»qeigen) 

c purpose: 

c This subroutine sets up the areas in the complex qll -plane 
c to search for the zeroes of the modal function, 
c 

c inputs... 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c outputs... 

c qeigen - complex array containing the zeros of the modal 
c function 

c nrmode - actual number of modes found 
c 

c subroutines called... 
c fzerox 

c findfx 

c common block areas... 
c coml 



mxlayr - maximum number of layers allowed in refractivity 
profile 

nzlayr - actual number of layers in refractivity profile 
aloss - maximum attenuation rate (in db/km) of modes 
to be found 

dm min - minimum of zim(j)-zim(l) 
tl - koal23 

dmesh - initial mesh size divisor 



implicit double precision (a-h,o-z) 
complex* 16 ctemp,tl,qeigen^eros,fxl,fx2 
parameter(c201og=8.68588963806504d0,one=( 1 .d0,0.d0), 
$ step0=1024.d0»step0h=step0/2.d0) 

character*3 filemjdine 
character*40 rfile 



qeigen 



c 
c 
c 
c 
c 
c 
c= 
c 
c 
c 
c 



- complex array containing all the zeros of the 
modal function found 

- complex array containing the zeros of the modal 
function found in the current rectangular region 
of the complex ql 1-plane 






42 

43 

44 

45 

46 

47 

48 Sinclude: ’mlaparm.inc’ 

***** Begin listing of: mlaparm.inc 



use include file for parameters of 
mxlayr max # layers 
mxmode max # modes 
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1 c 

2 c include file to define the 

3 c maximum # of layers (mxlayr) 

4 c maximum # of modes (mxmode) 

5 c 



6 parameter (mxlayr=35 ) 

7 parameter (mxmode=127 ) 
***** End listing of: mlaparm.inc 



49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 



dimension qeigen(mxmode),zeros(2*mxmode+l) 

common /coml/waveno 
c 

^itt^^^^^aiiait^^^tJti********************************************************** 

c rl - value of mode search on the real part of ql 1 at the 
c left edge in tmesh units. 

c r2 - value of mode search on the real part of qll at the 
c right edge in tmesh units, 

c bot - value of the imaginary part of qll at the bottom 
c edge (this is set to 0). 

c tO - value of the imaginary part of ql 1 at the top edge 
c in tmesh units, 

c step - size of search areas. 

c 

c 

c start of executable statements 

c — 

4tilt7lt:ti*4tHc**:tt:tt^t************************************ ******* ******** 

c set up search areas for finding modes and solve for modes, 
c calculate approximate value for re(qll) where modes start 
c occurring 

c 

nrmode=0 

c 

recons=dreal(tl) 
rstart=-2.0d-6*dm min 
qtest=rstart*recons 

ttopl=(2.0d-3/(waveno*c201og))*recons 
ttop=aloss*ttopl 
tmesh=ttop 1/dmesh 
if(tmesh .gt. O.ldO) tmesh=1.0d>l 
if(tmesh .gt. l.Od-3) then 
tol=1.0d-4 
else 

tol=tmesh*0.1d0 
end if 
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91 c 



92 c***** 


93 


write(*,1002) 


94 


write(16,1002) 


95 c 




96 


write(*,1003) tmesh.tol 


97 


write(16,1003) tmesh,toI 


98 c 




99 


tO=dnint(ttop/tmesh)+ 1 .dO 


100 


bot^O.dO 


101 


102 c 




103 


rl=dnint(qtest/tmesh)-stepOh 


104 


rl=rl 


105 


rr=rl+stepO 


106 


call rmdfx(rl,tO,fxl,xl,yl,lmesh) 


107 


il=int(rl) 


108 


in=int(stepO)+int(r 1) 


109 


nline=0 


110 50 


continue 


111 


do 100 nn=il,in 


112 


r2=rl+l.d0 


113 


call findfx(i2,t0,fx2,x2,y2,tmesh) 


114 


if (yl*y2 .It. O.dO) then 


115 


r=rl*tmesh 


116 


write (*.5555) r 


117 c 


write (16,5555) r 


118 5555 format (/5 x/f= ’, d28.16) 


119 


go to 400 


120 


endif 


121 


rl=r2 


122 


fxl=fx2 


123 


xl=x2 


124 


yl=y2 


125 100 


continue 


126 


return 


127 c 




128 400 


continue 


129 


niold=nrmode 


130 


nline=nline+l 


131 


nc 100=int(nline/100) 


132 


nc 10=int((nUne-nc 100* 100)/10) 


133 


nc l=nline-nc 100* 100-nc 10* 10 


134 


kline=char(48+nc 100)//char(48+nc 10)//char(48+nc 1) 


135 


rfile=’e:\ins5Nr7/rilem//kline//\dat’ 


136 c****** 


137 


open (unit=25,file=rfile,status=’unknown’) 


138 




139 


call fzerox(t0j 1 /1 1 ,fx 1 ,x 1 ,y 1 ,fx2,x2.y2,zeros,nrold,nmew, 


140 


$ tmesh) 
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141 close (25) 

142 c****** 

143 c 

144 nnnode=nmew 

145 newO=nrmode-nrold 

146 c 

147 Q********************’(‘****’*‘***********’*‘*********^**’^****************’*‘ 

148 c mode search completed 

149 c order zeros found by order of increasing real part. 

150 Q******************=<‘***********’*‘****’*‘’*‘****’*"*‘’^****’*"*‘’*‘*****’*'* ********’*'* 

151 c 

152 if(nrmode .gt, 1) then 

153 jkend=nrmode-l 

154 

155 do420jk=ljkend 

156 nend=nnnode-jk 

157 

158 do410 ja=l,nend 

159 nij=nrmode-ja 

160 mjl=nrj+l 

161 if(dimag(zeros(nrjl)).lt. 

162 $ dimag(zeros(nrj))) then 

1 63 ctemp=zeros(nij 1 ) 

164 zeros(nijl)=zeros(nij) 

165 zeros(nrj)=ctemp 

166 end if 

167 410 continue 

168 420 continue 

169 end if 

170 c 

171 Q********************************************************************** 

172 c the possibility exists that duplicate (within the tolerance ’tol’) 

173 c zeros of the modal equation will be found, eliminate these 

174 c duplicate zeros. 

175 Q************************************************ ********************** 

176 c 

177 jkflag=0 

178 jkend=nrmode-l 

179 

180 do240jk=ljkend 

181 jkl=jk+l 

182 ctemp=zerosOTcl) 

183 chksq=cdabs(zeros(jk-jkflag)-ctemp) 

184 if(chksq .It. tol) then 

185 jkflag=jkflag+l 

186 go to 240 

187 end if 

188 zerosO'kl-jkflag)=ctemp 

189 240 continue 

190 nrmode=nrmode-jkfIag 
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191 c 

193 c Store the zeros as the eigenvalues. 

194 Q*************************’*‘***************** 

195 c 

196 nrmode=minO(nrmode,mxmode) 

197 

198 do 430 jk= Ijirmode 

199 qeigen01c)=zeros(jk) 

200 430 continue 

201 c 

202 il=int(rll)+l 

203 if (il .It. in) then 

204 rl=rll+l.d0 

205 call findfx(rl,tO,fxKxUyl,tmesh) 

206 go to 50 

207 endif 

208 return 

209 c 

211 c format statements 

213 1000 format(/5x/searching for zeroes in this areas are’, 

214 $ ’ defined by:75x/istart= ’410/5x,’itop= \il0, 

215 $ 5x,’ ibot= ’411/5x,’ileft= \il0,5x/iright= \il0) 

216 

217 1001 format(5x44,’ new zeroes found in this area.’/) 

218 

219 1002 format(//5x/******* start search for modal eigenvalues’, 

220 $ ’ *♦***♦*’) 

221 



222 1003 format(/5x,’tmesh= ’,dl5.54x,’tol= ’,dl5.5/) 

223 

224 end 

225 c 

226 c 

227 

228 

229 c 

230 subroutine fzerox(t04'l,rll,fxl,xl,yli'x2,x2,y245eros4irold, 

231 $ nmew,tmsh0) 

232 c 



4i4t4t4i4i4i4c4< 

4< 4< 4< 4> 4> 4( 4( * 4< 4< 4< 4< 4c 4t 4> 4< 4> 4< 4< 4< 4< 4> 4> 4( 4( 4< 4< 4< 4< 4< 4< 4> 4> 4> 4t 4< 4< 4< 4< 4< 4< * 4< 4< 4< 4c 4t )ti 4c 4( 4( 4( 4i 



233 c***** 

234 c fzerox is a routine for finding the zeroes of a complex function, f, 

235 c which lie within a specified rectangular region of the 

236 c complex qll plane, assuming that the function has only 

237 c simple zeroes over this rectangle. 

238 c 

239 c parameters specifying the search rectangle: 

240 c tmesh - set equal to about half the average spacing between 
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241 c zeroes within the rectangle. A smaller value may be used 

242 c as a safety measure, but too small a value will result 

243 c in excessively long run time. 

244 c zeros - output list of (complex) values of ql 1 at which 

245 c zeroes are found. 

246 c nmew.nrold- the number of zeroes found 

247 c 

248 c subroutines called- 

249 c findfx 

250 c roots 

251 c***** 

252 implicit double precision (a-h,o-z) 

253 complex*16 fl0,f01,fll,fxl,fx2,fx00/xl0,fx01,fxll, 

254 + one,sol,zeros 

255 parameter(one=(l.d0,0.d0)) 

256 Sinclude: ’mlaparm.inc’ 

Begin listing of: mlaparm.inc 

1 c 

2 c include file to define the 

3 c maximum # of layers (mxlayr) 

4 c maximum # of modes (mxmode) 

5 c 

6 parameter (mxlayr=35 ) 

7 parameter (mxmode=127 ) 

♦♦♦♦♦ listing of: mlaparm.inc 

257 dimension sol(3),theta(2),zeros(2*mxmode+l),rline(2,1024) 

258 c 

259 common/tmccom/tmesh 

260 c***** 

261 c 

262 npointssO 

263 nmew=nrold 

264 tmesh=tmshO 

265 bot=0.d0 

266 t=t0-l.d0 

267 rll=rl 

268 r=rll 

269 c***** 

270 c 

271 fx01=fxl 

272 fxOlr^xl 

273 fx01i=yl 

274 fxll=0t2 

275 fxlli^x2 

276 fxlli=y2 

277 go to 82 

278 c 

279 c***** 

280 c enter mesh square from left side or exit rectangle at right edge. 

281 c 
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282 20 i^r+l.d0 

283 21 fx01=fxll 

284 fxOlr^fxllr 

285 fx01i=fxlli 

286 fxOO=fxlO 

287 fxOOr^fxlOr 

288 fxOOi=fxlOi 

289 22 continue 

290 call findfx(r+l.dO,t+l.dO,fxll,fxllrJxlli,tmesh) 

291 call findfx(r+l.dO,t^xlO/xlOr/xlOi,Unesh) 

292 c******* 

293 c Determine the edge of exit of im(f)=0 from current mesh. 

294 edgeit=fx01i*fxlli 

295 edgeib=fxOOi*fxlOi 

296 if (edgeib .gt. O.dO) then 

297 c lm(f)=0 goes through the 01 to 10 line. 

298 if (edgeit .gt. O.dO) then 

299 c lm(f)=0 goes through the 10 to 11 edge (edge 1). 

300 lout=l 

301 else 

302 c lm(0=0 goes through the 01 to 11 edge (edge 2) 

303 lout=2 

304 end if 

305 else 

306 c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

307 lout=4 

308 if (edgeit .iL O.dO) then 

309 c lm(0=0 also runs through 01 to 11 and 10 to 11 edges. 

310 c Store crossing location and in/out information. 

3 1 1 knot34=knot34+ 1 

312 c loc34r(knot34)=lr 

313 c loc34i(knot34)=li 

314 end if 

315 end if 

316 c******* 

317 go to 85 

318 c***** 

319 c enter mesh square from bottom side or exit rectangle at top edge. 

320 40 t=t+l.d0 

321 41 fx00=fx01 

322 fxOOi^fxOlr 

323 fx00i=fx01i 

324 fxl0=fxll 

325 fxlOi^fxllr 

326 fxlOUfxlli 

327 42 continue 

328 call findfx(r,t+l.d0»fx01^x01r,fx01i,tmesh) 

329 call fmdfx(r+ 1 .d0,t+ 1 .d0,fx 1 1 ,fx 1 1 r/x 1 1 i,tmesh) 

330 c******* 

331 c Determine the edge of exit of im(0=0 from current mesh. 
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332 edgeil=fx00i*fx01i 

333 edgeir^fxlOi*fxlli 

334 if (edgeir .gt. O.dO) then 

335 c lm(0=0 goes through the 00 to 1 1 line. 

336 if (edged .gt. 0.d0) then 

337 c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

338 lout=2 

339 else 

340 c lm(0=0 goes through the 00 to 01 edge (edge 3). 

341 lout=3 

342 end if 

343 else 

344 c lm(0=0 goes through the 10 to 1 1 edge (edge 1) 

345 lout=l 

346 if (edged .It. O.dO) then 

347 c lm(0=0 also runs through 00 to 01 and 01 to 11 edges. 

348 c Store crossing location and in/out information. 

349 knot4 l=knot4 1 + 1 

350 c loc41r(knot41)=lr 

351 c loc41i(knot41)=li 

352 end if 

353 end if 

354 c******* 

355 go to 85 

356 c***** 

357 c enter mesh square from right side or exit rectangle at left edge. 

358 

359 60 r=r-l.d0 

360 61 fxlUfxOl 

361 fxlln=fx01r 

362 fxlli=fx01i 

363 fxl0=fx00 

364 fxl0r=fx00r 

365 fxl0i=fx00i 

366 65 continue 

367 caU fmdfx(r.t+ 1 .dO/xO 1 JxO lr,fx0 1 i,tmesh) 

368 call findfx(r.ufx00,fx00r/x00i,tmesh) 

369 c******* 

370 c Determine the edge of exit of im(0=0 from current mesh. 

371 edgeit=fx01i*fxlli 

372 edgeib=fx00i*fxl0i 

373 if (edgeit .gt. O.dO) then 

374 c lm(0=0 goes through the 01 to 10 line. 

375 if (edgeib .gt. O.dO) then 

376 c lm(0=0 goes through the 00 to 01 edge (edge 3). 

377 lout=3 

378 else 

379 c lm(0=0 goes through the 00 to 10 edge (edge 4) 

380 lout=4 

381 end if 
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382 else 

383 c lm(0=0 goes through the 01 to 11 edge (edge 2) 

384 lout=2 

385 if (edgeib .It. O.dO) then 

386 c lm(0=0 also runs through 00 to 10 and 00 to 01 edges. 

387 c Store crossing location and in/out information. 

388 knot 12=knot 12+1 

389 c locl2r(knotl2)=lr 

390 c locl2i(knotl2)=li 

391 end if 

392 end if 

393 Q******* 

394 go to 85 

395 c***** 

396 c enter mesh square top side or exit rectangle at bottom edge. 

397 c 

398 80 t=t-l.d0 

399 81 fx01=fx00 

400 fxOlr^fxOOr 

401 fx01i=fx00i 

402 fxll=fxl0 

403 fxllr=fxl0r 

404 fxlli=fxl0i 

405 82 continue 

406 call findfx(r,t,fx00,fx00r/x00i,tmesh) 

407 call findfx(r+l.d0,t.fxl0jxl0r.fxl0i,tmesh) 

408 c 

409 Q******* 

410 c Determine the edge of exit of im(0=0 from current mesh. 

411 c 

412 83 continue 

413 edgeil=fx00i*fx01i 

414 edgeir=fxl0i*fxlli 

415 if (edged .gt. O.dO) then 

416 c lm(0=0 goes through the 00 to 11 line. 

417 if (edgeir .gt. O.dO) then 

418 c lm(0=0 goes through the 00 to 10 edge (edge 4) 

419 lout=4 

420 else 

421 c lm(0=0 goes through the 10 to 11 edge (edge 1). 

422 lout=l 

423 end if 

424 else 

425 c lm(f)=0 goes through the 00 to 01 edge (edge 3) 

426 lout=3 

427 if (edgeir .It. O.dO) then 

428 c lm(0=0 also runs through 00 to 10 and 10 to 11 edges. 

429 c Store crossing location and in/out information. 

430 knot23=knot23+l 

431 c loc23r(knot23)=lr 
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432 c loc23i(knot23)=li 

433 end if 

434 end if 

435 c*************** *************************************************** ****** 

436 c Test for there being at least one re(f)=0 line entering and 

437 c leaving the mesh square. 

438 Q*’*'**’*'***********’^*********************************************** ******** 

439 c 

440 85 continue 

441 npoint=npoint+l 

442 nr=r*tmesh 

443 ttt=t*tmesh 

444 rline( 1 4 ipoint)=rrr 

445 rline( 24 ipoint)=ttt 

446 if (npoint ,eq. 1024) lout=5 

447 c 

448 if((t .It. bot) .or, (t .gt. tO)) go to 100 

449 

450 if ((fx00r*fxl0r .gt. O.dO) .and. (fx01r*fxl Ir ,gt. O.dO) 

451 + .and, (fx00r*fx01r .gt. O.dO)) go to (20,40,60,80,100) 

452 + lout 

453 Q**********************^************************************************** 

454 c Computate the values of the modal function at the comers of 

455 c a mesh square to determine its Taylor series to the 3rd order 

456 c for estimating its root locations, 

457 Q********************************************************************** 

458 c 

459 c f00=one 

460 fl0=cdexp(fxl0*fx00)-one 

461 f01=cdexp(fx01-fx00)-one 

462 fl l=cdexp(fxl l-fx00)-one 

463 c 

464 c*********** estimate locations of zeroes by radicals***************** 

465 c 

466 call roots(f 10401 41 1 ,sol jirsol) 

467 c 

468 do 90 n=l,nrsol 

469 ureal = dreal(sol(n)) 

470 uimag = dimag(sol(n)) 

471 if (ureal .It. O.dO .or. ureal .gt. l.OdO) go to 90 

472 if (uimag .It, O.dO .or. uimag .gt, l.OdO) go to 90 

473 92 theta(l)=(r+ureal)*tmesh 

474 theta(2)=(t+uimag)*tmesh 

475 nmew = nmew+1 

476 zeros(nmew)=dcmpbc(theta( l),theta(2)) 

477 90 continue 

478 c 

479 

480 c continue following the phase line 

481 c******************************************* 
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482 c 

483 95 continue 

484 go to (20,40,60,80,100) lout 

485 c** 

486 100 continue 

487 c****** 

488 write (25,8888) ((rline(i,j), i=l, 2), j=l, npoint) 

489 8888 format (2e 15.5) 

490 c 

491 return 

492 end 
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APPENDIX C: CONSTANT PHASE LINE TRACKING FLOW CHARTS 
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APPENDIX D: CONSTANT PHASE LINES IN THE COMPLEX PLANE 
This appendix contains plots of the constant phase lines Im{D(q)} = 0 and 
Re{D(q)} = 0 initiating from the top edge of the search region in the complex q^^/tj 
plane for the evaporation ducts of 2, 4, 6, 8, 10, 20, 30, and 40 m heights. 
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Figure D.l Constant phase lines in the qu/tj plane (2 m duct). 
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Figure D.3 Constant phase lines in the plane (6 m duct). 
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xlO-6 Constant phase lines in the ql 1/tl plane for 10 in duct 
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xlO'6 Constant phase lines in the ql 1/tl plane for 30 m duct 




Figure D.7 Constant phase lines in the q^j/tj plane (30 m duct). 
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Figure D.8 Constant phase lines in the qjj/tj plane (40 m duct). 
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APPENDIX E: SEPARATION OF Im{D(q))=0 LINES 



This appendix contains figures of the separation of Im{D(q)}=0 lines along the top 
edge of the search region. The 2, 4, 6, 8, 10, 20, 30 and 40 meter evaporation ducts are 
included. 
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Figure E.1 Separation of Im{D(q))=0 lines along the top edge of the search region 
(2 m duct). 
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Figure E.2 Separation of lm{D(q))=0 lines along the top edge of the search region 
(4 m duct). 
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Figure E.3 Separation of Im{D(q))=0 lines along the top edge of the search region 
(6 m duct). 
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(8 m duct). 
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Rgurc E.5 Separation of Im{D(q))=0 lines along the top edge of the search region 
(10 m duct). 
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(20 m duct). 
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Figure E.7 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(30 m duct). 
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Figure E.8 Separation of Im{D(q))=0 lines along the top edge of the search region 
(40 m duct). 
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APPENDIX F: DIFFERENCE BETWEEN CONSECUTIVE Re(qytj) VALUES 
Separation in real part of neighboring q-eigenvalues scaled by t^ is plotted in this 
appendix against the eigenvalue locations along the real q^ ^/t ^ axis. 
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Figure F.l Difference between consecutive ReCq^^tj) values (2 m duct). 
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Figure F.2 Difference between consecutive Re(qj^|) values (4 m duct). 
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Figure F.3 Difference between consecutive values (6 m duct). 
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Figure F.4 Difference between consecutive R^Cq^^tj) values (8 m duct). 
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Figure F.5 Difference between consecutive values (10 m duct). 
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Figure F.6 Difference between consecutive R^Cq^^tj) values (20 m duct). 
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Figure F.7 Difference between consecutive Re(qj^tj) values (30 m duct). 
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Figure F.8 Difference between consecutive Re(qj^tj) values (40 m duct). 
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